#include<cstring>
#include<time.h>
using namespace std;
#include "../code/linear_algebra.h"
#include "../code/HyperAlgebra.h"
#include "../code/random.h"
//#include"../../HyperLex/extern_head/HyperLex.h"
#include"../extern/header/lexical.h"
int test_entrance(const char* output_path);
static int cores_amount_from_input(const char* s);
int main(int argc, char* argv[])
{
	int info;
	//parameter_file pf;
	int amount_of_cores;
	const char* output_path;
	output_path = argv[1];
	amount_of_cores = cores_amount_from_input(argv[2]);
	info = 0;
	info = test_entrance(output_path);
	return info;
}







static int cores_amount_from_input(const char* s)
{
	hyperlex::CharBuffer temp;
	int core;
	temp << s;
	core = temp.DequeueInt();
	return core;
}
static int test_01(const char* output_path, hyperlex::ParaFile& pf);
static int test_02(const char* output_path, hyperlex::ParaFile& pf);
static int test_03(const char* output_path, hyperlex::ParaFile& pf);
static int test_04(const char* output_path, hyperlex::ParaFile& pf);
int test_entrance(const char* output_path)
{
	int item;
	int info;
	FILE* fp;
	size_t site;
	hyperlex::ParaFile pf;
	pf.initial("./parameter/test_item.txt");
	//===========================================================
	site = pf.SearchKey("item", hyperlex::ParaFile::Int);
	if (site != pf.Amount())
		item = pf.FirstInt(site);
	else 
		fprintf(stderr, "No item!\n");
	cout << "item: " << item << endl;
	info = 0;
	fp = stdout;
	fprintf(stdout, "============================================\n");
	switch (item)
	{
	case 1:
		fprintf(stdout, "Test inner dot, and Norm2\n");
		info = test_01(output_path, pf);
		break;
	case 2:
		fprintf(stdout, "Test gemv.\n");
		info = test_02(output_path, pf);
		break;
	case 3:
		fprintf(stdout, "Test Substract, Norm2Normalization, FiniteDifference.\n");
		info = test_03(output_path, pf);
		break;
	case 4:
		info = test_04(output_path, pf);
		break;
	default:
        info = test_01(output_path, pf);
		break;
	}
	fprintf(stdout, "============================================\n");
	return info;
}
static void test_inner(size_t length);
static int test_01(const char* output_path, hyperlex::ParaFile& pf)
{
	double A[3], B[3];
	HyperAlgebra::AlgebraMKL mkl;
	HyperAlgebra::BasicVector BV;
	
	A[0] = 0.1;
	A[1] = 1.0;
	A[2] = 0.0;
	B[0] = 1.0;
	B[1] = 0.1;
	B[2] = 1.0;
	printf("BV:%lf\n", BV.InnerDot(A, B, 3));
	printf("mkl:%lf\n", mkl.InnerDot(A, B, 3));
	test_inner(3);
	test_inner(8);
	test_inner(32);
	test_inner(122);
	test_inner(256);
	test_inner(1023);
	test_inner(1024*32);
}
static void test_gemv(size_t row, size_t column, size_t deltaLd);
static int test_02(const char* output_path, hyperlex::ParaFile& pf)
{
	test_gemv(1, 1, 0);
	test_gemv(1, 1, 1);
	test_gemv(4, 1, 0);
	test_gemv(4, 1, 1);
	test_gemv(1, 4, 0);
	test_gemv(1, 4, 1);

	test_gemv(32, 32, 0);
	test_gemv(32, 32, 1);
	test_gemv(256, 256, 0);
	test_gemv(256, 256, 1);
	test_gemv(1, 256, 0);
	test_gemv(1, 256, 1);
	test_gemv(256, 1, 0);
	test_gemv(256, 1, 1);
}
static void test_vec2(size_t length);
static int test_03(const char* output_path, hyperlex::ParaFile& pf)
{
	test_vec2(1);
	test_vec2(2);
	test_vec2(5);
	test_vec2(8);
	test_vec2(15);
	test_vec2(32);
	test_vec2(256);
	test_vec2(256 * 256);
}

static void test_sysv(size_t dim, size_t deltaLd);
static int test_04(const char* output_path, hyperlex::ParaFile& pf)
{
	double A[9];
	double W[3];
	size_t i, j;
	HyperAlgebra::AlgebraMKL mkl;

	A[0] = 4; A[1] = 2; A[2] = 2;
	A[3] = 2; A[4] = 4; A[5] = 2;
	A[6] = 2; A[7] = 2; A[8] = 4;
	for (i = 0; i < 3; i++)
	{
		for (j = 0; j < 3; j++)
		{
			printf("%lf ", A[i * 3 + j]);
		}
		printf("\n");
	}
	printf("===================================\n");
	mkl.MKLdsyev(HyperAlgebra::lower, 3, A, 3, W);
	for (i = 0; i < 3; i++)
	{
		for (j = 0; j < 3; j++)
		{
			printf("%lf ", A[i * 3 + j]);
		}
		printf("\n");
	}
	printf("===================================\n");
	for (i = 0; i < 3; i++) printf("%lf ", W[i]);
	test_sysv(1, 0);
	test_sysv(1, 1);
	test_sysv(2, 0);
	test_sysv(2, 1);
	test_sysv(15, 0);
	test_sysv(15, 1);
	test_sysv(256, 0);
	test_sysv(256, 1);
	test_sysv(256, 257);
	test_sysv(256 * 16, 0);
	test_sysv(256 * 16, 1);
	test_sysv(256 * 16, 257);
}



static void test_inner(size_t length)
{
	HyperAlgebra::AlgebraMKL mkl;
	HyperAlgebra::BasicVector BV;
	HYPER_random Hr;
	double *A, *B;
	double a, b, c, d;
	size_t i;
	A = (double*)malloc(sizeof(double) * length);
	B = (double*)malloc(sizeof(double) * length);
	for (i = 0; i < length; i++)
	{
		A[i] = Hr.rand(-1.0, 1.0);
		B[i] = Hr.rand(-1.0, 1.0);
	}
	a = BV.InnerDot(A, B, length);
	b = mkl.InnerDot(A, B, length);
	printf("length: %llu\n", length);
	printf("\tBV.InnerDot(A, B, length):%lf, mkl.InnerDot(A, B, length):%lf delta: %.16lf\n", a, b, a - b);
	a = BV.InnerDot(A, A, length);
	b = mkl.InnerDot(A, A, length);
	c = BV.Norm2(A, length);
	d = mkl.Norm2(A, length);
	printf("\tBV.InnerDot(A, A, length):%lf, mkl.InnerDot(A, A, length):%lf delta: %.16lf\n", a, b, a - b);
	printf("\tBV.Norm2(A, length):%lf, mkl.Norm2(A, length):%lf delta: %.16lf\n", c, d, c - d);

	free(A);
	free(B);
}
static void test_gemv(size_t row, size_t column, size_t deltaLd)
{
	HyperAlgebra::AlgebraMKL mkl;
	HyperAlgebra::BasicMatrix BM;
	HyperAlgebra::BasicVector BV;
	HYPER_random Hr;
	double *A, *X;
	double *Y_1, *Y_2;
	size_t i;
	size_t ldM;
	double a, b, c;
	ldM = column + deltaLd;
	printf("row: %llu, column: %llu, ldM: %llu\n", row, column, ldM);
	A = (double*)malloc(sizeof(double) * row * ldM);
	X = (double*)malloc(sizeof(double) * column);
	Y_1 = (double*)malloc(sizeof(double) * row);
	Y_2 = (double*)malloc(sizeof(double) * row);
	for(i = 0; i < row * ldM; i++) A[i] = Hr.rand(-1.0, 1.0);
	for(i = 0; i < column; i++) X[i] = Hr.rand(-1.0, 1.0);
	BM.GEMVn(A, ldM, X, Y_1, row, column);
	mkl.GEMVn(A, ldM, X, Y_2, row, column);

	a = BV.Norm2(Y_1, row);
	b = BV.Norm2(Y_2, row);
	c = BV.Norm2(Y_1, Y_2, row);
	printf("\tBM.GEMVn: %lf, mkl.GEMVn: %lf, BM.GEMVn - mkl.GEMVn: %.16lf\n", a, b, c);

	free(X);
	free(Y_1);
	free(Y_2);
	X = (double*)malloc(sizeof(double) * row);
	Y_1 = (double*)malloc(sizeof(double) * column);
	Y_2 = (double*)malloc(sizeof(double) * column);
	for(i = 0; i < row; i++) X[i] = Hr.rand(-1.0, 1.0);
	BM.GEMVt(A, ldM, X, Y_1, row, column);
	mkl.GEMVt(A, ldM, X, Y_2, row, column);
	
	a = BV.Norm2(Y_1, column);
	b = BV.Norm2(Y_2, column);
	c = BV.Norm2(Y_1, Y_2, column);
	printf("\tBM.GEMVn: %lf, mkl.GEMVn: %lf, BM.GEMVn - mkl.GEMVn: %.16lf\n", a, b, c);


	free(A);
	free(X);
	free(Y_1);
	free(Y_2);



}
static void test_vec2(size_t length)
{
	HyperAlgebra::AlgebraMKL mkl;
	HyperAlgebra::BasicVector BV;
	HYPER_random Hr;
	double *vector;
	double *vector_1, *vector_2;
	double *Y_1, *Y_2;
	size_t i;
	double a;
	double DeltaReverse;
	printf("length: %llu\n", length);
	vector = (double*)malloc(sizeof(double) * length);
	vector_1 = (double*)malloc(sizeof(double) * length);
	vector_2 = (double*)malloc(sizeof(double) * length);
	Y_1 = (double*)malloc(sizeof(double) * length);
	Y_2 = (double*)malloc(sizeof(double) * length);
	for(i = 0; i < length; i++) 
	{
		vector[i] = Hr.rand(-1.0, 1.0);
		vector_2[i] = Hr.rand(-1.0, 1.0);
		Y_1[i] = vector[i];
		Y_2[i] = vector[i];
	}
	if(length < 16)
	{
		printf("\tvector: ");
		for(i = 0; i < length; i++)  printf("%lf  ", vector[i]);
		printf("\n");
	}
	printf("\tNorm2 input: %.16lf\n", BV.Norm2(vector, length));
	BV.Norm2Normalization(Y_1, length);
	mkl.Norm2Normalization(Y_2, length);

	printf("\tBV.Norm2: %.16lf, mkl.Norm2: %.16lf", BV.Norm2(Y_1, length), BV.Norm2(Y_2, length));
	printf(", a - b: %.16lf\n", BV.Norm2(Y_1, Y_2, length));
	
	if(length < 16)
	{
		printf("\tY_1: ");
		for(i = 0; i < length; i++)  printf("%lf  ", Y_1[i]);
		printf("\n");
	}
	//====================================
	for(i = 0; i < length; i++) 
	{
		Y_1[i] = vector[i];
		Y_2[i] = vector[i];
	}
	a = BV.Norm2(vector_2, length);
	a = a * a;
	BV.Substract(Y_1, length, vector_2, a);
	mkl.Substract(Y_2, length, vector_2, a);

	printf("\tBV.Substract: %.16lf, mkl.Substract: %.16lf", BV.Norm2(Y_1, length), BV.Norm2(Y_2, length));
	printf(", a - b: %.16lf\n", BV.Norm2(Y_1, Y_2, length));

	if(length < 16)
	{
		printf("\tY_1: ");
		for(i = 0; i < length; i++)  printf("%lf  ", Y_1[i]);
		printf("\n");
	}
	//====================================
	DeltaReverse = 1.0E-6;

	for(i = 0; i < length; i++)
	{
		vector_1[i] = vector[i] + DeltaReverse * vector_2[i];
	}
	BV.FiniteDifference(Y_1, length, 1.0 / DeltaReverse, vector_1, vector);
	mkl.FiniteDifference(Y_2, length, 1.0 / DeltaReverse, vector_1, vector);
	printf("\tBV.Substract: %.16lf, mkl.Substract: %.16lf", BV.Norm2(Y_1, length), BV.Norm2(Y_2, length));
	printf(", a - b: %.16lf\n", BV.Norm2(Y_1, Y_2, length));
	printf("\tE BV.Substract: %.16lf, E mkl.Substract: %.16lf\n", BV.Norm2(Y_1, vector_2, length), BV.Norm2(Y_2, vector_2, length));

	free(vector);
	free(vector_1);
	free(vector_2);
	free(Y_1);
	free(Y_2);
}
static void test_sysv(size_t dim, size_t deltaLd)
{
	HyperAlgebra::AlgebraMKL mkl;
	HyperAlgebra::BasicMatrix BM;
	HyperAlgebra::BasicVector BV;
	HYPER_random Hr;
	double *A, *eigenvalues;
	double *H, *temp, *Y;
	size_t i, j;
	size_t ldM;
	double a, b, c;
	ldM = dim + deltaLd;
	printf("row: %llu, column: %llu, ldM: %llu\n", dim, dim, ldM);
	A = (double*)malloc(sizeof(double) * dim * ldM);
	eigenvalues = (double*)malloc(sizeof(double) * dim);
	H = (double*)malloc(sizeof(double) * dim * ldM);
	temp = (double*)malloc(sizeof(double) * dim * dim);
	Y = (double*)malloc(sizeof(double) * dim * dim);
	for(i = 0; i < dim; i++)
	{
		A[i * ldM + i] = Hr.rand(-1.0, 1.0);
		for (j = i + 1; j < dim; j++)
		{
			A[i * ldM + j] = Hr.rand(-1.0, 1.0);
			A[j * ldM + i] = A[i * ldM + j];
		}
	} 
		
	for(i = 0; i < dim; i++)
		for (j = 0; j < dim; j++) H[i * ldM + j] = A[i * ldM + j];
	mkl.MKLdsyev(HyperAlgebra::upper, dim, H, ldM, eigenvalues);
	for(i = 0; i < dim * dim; i++) temp[i] = 0.0;
	for (i = 0; i < dim; i++)
		temp[i * dim + i] = eigenvalues[i];
	mkl.gemmNN(dim, dim, dim, 1.0, H, ldM, temp, dim, 0.0, Y, dim);
	mkl.gemmNT(dim, dim, dim, 1.0, Y, dim, H, ldM, 0.0, temp, dim);
	printf("\tHyperAlgebra::upper Norm Max: %.16lf\n", BM.NormMax(dim, dim, temp, dim, A, ldM));

	for(i = 0; i < dim; i++)
		for (j = 0; j < dim; j++) H[i * ldM + j] = A[i * ldM + j];
	mkl.MKLdsyev(HyperAlgebra::lower, dim, H, ldM, eigenvalues);
	for(i = 0; i < dim * dim; i++) temp[i] = 0.0;
	for (i = 0; i < dim; i++)
		temp[i * dim + i] = eigenvalues[i];
	mkl.gemmNN(dim, dim, dim, 1.0, H, ldM, temp, dim, 0.0, Y, dim);
	mkl.gemmNT(dim, dim, dim, 1.0, Y, dim, H, ldM, 0.0, temp, dim);
	printf("\tHyperAlgebra::lower Norm Max: %.16lf\n", BM.NormMax(dim, dim, temp, dim, A, ldM));

	free(A);
	free(eigenvalues);
	free(H);
	free(temp);
	free(Y);
}


